#include "number.h"

mpz_t tmp_complex_print_int;

mpfr_t tmp_complex_from_float_str,
       tmp_divide_re,
       tmp_divide_im,
       tmp_abs_n_re,
       tmp_abs_ll_re,
       tmp_abs_lu_re,
       tmp_abs_ul_re,
       tmp_abs_uu_re,
       tmp_abs_re_abs,
       tmp_abs_im_abs,
       tmp_abs_zero,
       tmp_diagonal,
       tmp_sqrt_n_re,
       tmp_sqrt_n_im;

mpc_t tmp_divide_numerator,
      tmp_divide_denominator,
      tmp_divide_quotient,
      tmp_abs_n,
      tmp_abs_ll,
      tmp_abs_lu,
      tmp_abs_ul,
      tmp_abs_uu,
      tmp_sqrt_n,
      tmp_sqrt_ll,
      tmp_sqrt_lu,
      tmp_sqrt_ul,
      tmp_sqrt_uu,
      tmp_sqrt_re_sqrt;

void init_complex_tmp_nums ()
{
    mpz_init(tmp_complex_print_int);
    mpfr_init2(tmp_complex_from_float_str, PREC);
    mpfr_init2(tmp_divide_re, PREC);
    mpfr_init2(tmp_divide_im, PREC);
    mpfr_init2(tmp_abs_ll_re, PREC);
    mpfr_init2(tmp_abs_n_re, PREC);
    mpfr_init2(tmp_abs_ll_re, PREC);
    mpfr_init2(tmp_abs_lu_re, PREC);
    mpfr_init2(tmp_abs_ul_re, PREC);
    mpfr_init2(tmp_abs_uu_re, PREC);
    mpfr_init2(tmp_abs_re_abs, PREC);
    mpfr_init2(tmp_abs_im_abs, PREC);
    mpfr_init2(tmp_abs_zero, PREC);
    mpfr_init2(tmp_diagonal, PREC);
    mpfr_init2(tmp_sqrt_n_re, PREC);
    mpfr_init2(tmp_sqrt_n_im, PREC);
    
    mpc_init2(tmp_divide_numerator, PREC);
    mpc_init2(tmp_divide_denominator, PREC);
    mpc_init2(tmp_divide_quotient, PREC);
    mpc_init2(tmp_abs_n, PREC);
    mpc_init2(tmp_abs_ll, PREC);
    mpc_init2(tmp_abs_lu, PREC);
    mpc_init2(tmp_abs_ul, PREC);
    mpc_init2(tmp_abs_uu, PREC);
    mpc_init2(tmp_sqrt_n, PREC);
    mpc_init2(tmp_sqrt_ll, PREC);
    mpc_init2(tmp_sqrt_lu, PREC);
    mpc_init2(tmp_sqrt_ul, PREC);
    mpc_init2(tmp_sqrt_uu, PREC);
    mpc_init2(tmp_sqrt_re_sqrt, PREC);
}

/*
 *  Sets the value of num to the complex value represented by re and im.
 */
void mpc_set_r_r (mpc_t num, Real* re, Real* im)
{
    char buf[10000];
    char tmp_re_buf[10000];
    char tmp_im_buf[10000];
    char re_buf[10000];
    char im_buf[10000];
    
    if (REAL(re))
    {
        gmp_sprintf(tmp_re_buf, "%+Zd", re->num);
        sprintf(re_buf, "%c0.%se%ld", *tmp_re_buf, tmp_re_buf + 1, re->exp);
    }
    else
    {
        sprintf(re_buf, "%s", NaN(re) ? "@nan@" : POS_INF(re) ? "@inf@" : "-@inf@");
    }
    
    if (!im)
    {
        sprintf(im_buf, "0");
    }
    else if (REAL(im))
    {
        gmp_sprintf(tmp_im_buf, "%+Zd", im->num);
        sprintf(im_buf, "%c0.%se%ld", *tmp_im_buf, tmp_im_buf + 1, im->exp);
    }
    else
    {
        sprintf(im_buf, "%s", NaN(im) ? "@nan@" : POS_INF(im) ? "@inf@" : "-@inf@");
    }
    
    sprintf(buf, "(%s %s)", re_buf, im_buf);
    mpc_set_str(num, buf, 10, MPC_RNDNN);
}

/*
 *  Sets the value of num to the floating point value represented by real.
 */
void mpfr_set_r (mpfr_t num, Real* real)
{
    char buf[10000];
    char real_buf[10000];
    
    if (REAL(real))
    {
        gmp_sprintf(real_buf, "%+Zd", real->num);
        sprintf(buf, "%c0.%se%ld", *real_buf, real_buf + 1, real->exp);
    }
    else
    {
        sprintf(buf, "%s", NaN(real) ? "@nan@" : POS_INF(real) ? "@inf@" : "-@inf@");
    }
    
    mpfr_set_str(num, buf, 10, MPFR_RNDN);
}

/*
 *  Frees the memory block pointed to by complex.
 */
void complex_free (Complex* complex)
{
    if (complex)
    {
        if (complex->re)
        {
            interval_free(complex->re);
        }
        
        if (complex->im)
        {
            interval_free(complex->im);
        }
        
        free(complex);
        complex = NULL;
    }
}

/*
 *  Returns the size of the string representation of complex as a long.
 */
long complex_strlen (Complex* complex)
{
    if (interval_zero_p(complex->im))
    {
        return interval_strlen(complex->re);
    }
    else
    {
        int sign = 1;
        int i = 1;
        int re = interval_strlen(complex->re);
        int im = interval_strlen(complex->im);
        
        return re + sign + im + i;
    }
}

/*
 *  Puts a string representation of complex into the string pointed to by buf.
 *  
 *    complex: [1:1 < (2:1) < 3:1]+[0:0 ≤ (0:0) ≤ 0:0]i => str: "[1:1 < (2:1) < 3:1]"
 *    complex: [1:1 < (2:1) < 3:1]+[1:0 ≤ (1:0) ≤ 1:0]i => str: "[1:1 < (2:1) < 3:1]+1:0i"
 */
void complex_print (char* buf, Complex* complex)
{
    if (interval_zero_p(complex->im))
    {
        interval_print(buf, complex->re);
    }
    else
    {
        char re_buf[interval_strlen(complex->re)];
        char im_buf[interval_strlen(complex->im)];
        
        interval_print(re_buf, complex->re);
        interval_print(im_buf, complex->im);
        
        sprintf(buf, "%s+%si", re_buf, im_buf);
    }
}

/*
 *  Puts a string representation of the real part of complex as a truncated
 *  integer into the string pointed to by buf.
 */
void complex_print_int (char* buf, Complex* complex)
{
    Real* rounded = real_round_origin(complex->re->n);
    
    if (!REAL(rounded))
    {
        sprintf(buf, "%s", NaN(complex->re->n) ? "NaN" : POS_INF(complex->re->n) ? "Infinity" : "-Infinity");
    }
    else if (ZERO(rounded))
    {
        sprintf(buf, "0");
    }
    else
    {
        mpz_set(tmp_complex_print_int, rounded->num);
        num_pad(tmp_complex_print_int, rounded->exp - num_len(tmp_complex_print_int));
        mpz_get_str(buf, 10, tmp_complex_print_int);
    }
    
    real_free(rounded);
}

/*
 *  Puts a string representation of the real part of complex as a fractional
 *  decimal into the string pointed to by buf.
 */
void complex_print_decimal (char* buf, Complex* complex)
{
    char real_buf[10000];
    
    if (REAL(complex->re->n))
    {
        gmp_sprintf(real_buf, "%+Zd", complex->re->n->num);
        sprintf(buf, "%c0.%se%ld", *real_buf, real_buf + 1, complex->re->n->exp);
    }
    else
    {
        sprintf(buf, "%s", NaN(complex->re->n) ? "NaN" : POS_INF(complex->re->n) ? "Infinity" : "-Infinity");
    }
}

/*
 *  Returns a pointer to a newly allocated Complex with the given value. The
 *  given intervals will be referenced by the new complex and must be freed only
 *  when this complex is freed.
 *  
 *    re: {1:1 < 2:1 < 3:1}, im: {0:1 ≤ 0:1 ≤ 0:1} => complex: {1:1 < 2:1 < 3:1}, {0:1 ≤ 0:1 ≤ 0:1}
 */
Complex* complex_new (Interval* re, Interval* im)
{
    Complex* complex = malloc(sizeof(Complex));
    
    if (complex == NULL)
    {
        return NULL;
    }
    
    memset(complex, 0, sizeof(Complex));
    complex->re = re;
    complex->im = im;
    
    return complex;
}

/*
 *  Returns a pointer to a newly allocated Complex with the given real as its
 *  real part and 0 as its imaginary part.
 */
Complex* complex_from_real (Real* real)
{
    Interval* re = interval_from_real(real);
    Interval* im = interval_zero();
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex with the integer represented
 *  by str as its real part.
 */
Complex* complex_from_int_str (char* buf)
{
    Real* real = real_from_str(buf, strlen(buf) - (*buf == '-'));
    return complex_from_real(real);
}

/*
 *  Returns a pointer to a newly allocated Complex with the float represented by
 *  str as its real part.
 */
Complex* complex_from_float_str (char* buf)
{
    mpfr_set_str(tmp_complex_from_float_str, buf, 10, MPFR_RNDN);
    Real* real = real_from_mpfr(tmp_complex_from_float_str, ROUND_NEAREST);
    
    return complex_from_real(real);
}

/*
 *  Returns a pointer to a newly allocated Complex containing a NaN value.
 */
Complex* complex_nan ()
{
    Interval* re = interval_nan();
    Interval* im = interval_zero();
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing only +∞.
 */
Complex* complex_pos_inf ()
{
    Interval* re = interval_pos_inf();
    Interval* im = interval_zero();
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing only -∞.
 */
Complex* complex_neg_inf ()
{
    Interval* re = interval_neg_inf();
    Interval* im = interval_zero();
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing only 0.
 */
Complex* complex_zero ()
{
    Interval* re = interval_zero();
    Interval* im = interval_zero();
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing only 1.
 */
Complex* complex_one ()
{
    Interval* re = interval_one();
    Interval* im = interval_zero();
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing π, accurate to
 *  digs digits.
 */
Complex* complex_pi ()
{
    Interval* re = interval_pi();
    Interval* im = interval_zero();
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing e, accurate to
 *  digs digits.
 */
Complex* complex_e ()
{
    Interval* re = interval_e();
    Interval* im = interval_zero();
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing only i.
 */
Complex* complex_i ()
{
    Interval* re = interval_zero();
    Interval* im = interval_one();
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex with the same value as the
 *  given complex.
 */
Complex* complex_dup (Complex* complex)
{
    Interval* re = interval_dup(complex->re);
    Interval* im = interval_dup(complex->im);
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing a + b.
 */
Complex* complex_add (Complex* a, Complex* b)
{
    Interval* re = interval_add(a->re, b->re);
    Interval* im = interval_add(a->im, b->im);
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing a - b.
 */
Complex* complex_subtract (Complex* a, Complex* b)
{
    Interval* re = interval_subtract(a->re, b->re);
    Interval* im = interval_subtract(a->im, b->im);
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing a * b.
 */
Complex* complex_multiply (Complex* a, Complex* b)
{
    Interval* re_re = interval_multiply(a->re, b->re);
    Interval* re_im = interval_multiply(a->re, b->im);
    Interval* im_re = interval_multiply(a->im, b->re);
    Interval* im_im = interval_multiply(a->im, b->im);
    
    Interval* re = interval_subtract(re_re, im_im);
    Interval* im = interval_add(re_im, im_re);
    
    interval_free(re_re);
    interval_free(re_im);
    interval_free(im_re);
    interval_free(im_im);
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing a / b.
 */
Complex* complex_divide (Complex* a, Complex* b)
{
    Interval* re;
    Interval* im;
    
    Interval* denominator_re = interval_multiply(b->re, b->re);
    Interval* denominator_im = interval_multiply(b->im, b->im);
    Interval* denominator = interval_add(denominator_re, denominator_im);
    
    if (real_zero_p(denominator->l) || real_zero_p(denominator->u) || interval_span_zero_p(denominator))
    {
        mpfr_set_prec(tmp_divide_re, PREC);
        mpfr_set_prec(tmp_divide_im, PREC);
        mpc_set_prec(tmp_divide_quotient, PREC);
        mpc_set_prec(tmp_divide_numerator, PREC);
        mpc_set_prec(tmp_divide_denominator, PREC);
        
        mpc_set_r_r(tmp_divide_numerator, a->re->n, a->im->n);
        mpc_set_r_r(tmp_divide_denominator, b->re->n, b->im->n);
        
        mpc_div(tmp_divide_quotient, tmp_divide_numerator, tmp_divide_denominator, MPC_RNDNN);
        
        mpc_real(tmp_divide_re, tmp_divide_quotient, MPFR_RNDN);
        mpc_imag(tmp_divide_im, tmp_divide_quotient, MPFR_RNDN);
        
        re = interval_new(real_neg_inf(), true, real_from_mpfr(tmp_divide_re, ROUND_NEAREST), true, real_pos_inf());
        im = interval_new(real_neg_inf(), true, real_from_mpfr(tmp_divide_im, ROUND_NEAREST), true, real_pos_inf());
    }
    else
    {
        Interval* re_re = interval_multiply(a->re, b->re);
        Interval* re_im = interval_multiply(a->re, b->im);
        Interval* im_re = interval_multiply(a->im, b->re);
        Interval* im_im = interval_multiply(a->im, b->im);
        
        Interval* numerator_re = interval_add(re_re, im_im);
        Interval* numerator_im = interval_subtract(im_re, re_im);
        
        re = interval_divide(numerator_re, denominator);
        im = interval_divide(numerator_im, denominator);
        
        interval_free(re_re);
        interval_free(re_im);
        interval_free(im_re);
        interval_free(im_im);
        interval_free(numerator_re);
        interval_free(numerator_im);
    }
    
    interval_free(denominator_re);
    interval_free(denominator_im);
    interval_free(denominator);
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing a % b.
 */
Complex* complex_modulo (Complex* a, Complex* b)
{
    Interval* re = interval_zero();
    Interval* im = interval_zero();
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing the absolute value
 *  of complex.
 */
Complex* complex_abs (Complex* complex)
{
    Real* l_abs;
    Real* u_abs;
    
    mpfr_set_prec(tmp_abs_n_re, PREC);
    mpfr_set_prec(tmp_abs_ll_re, PREC);
    mpfr_set_prec(tmp_abs_lu_re, PREC);
    mpfr_set_prec(tmp_abs_ul_re, PREC);
    mpfr_set_prec(tmp_abs_uu_re, PREC);
    mpfr_set_prec(tmp_abs_re_abs, PREC);
    mpfr_set_prec(tmp_abs_im_abs, PREC);
    mpfr_set_prec(tmp_abs_zero, PREC);
    mpc_set_prec(tmp_abs_n, PREC);
    mpc_set_prec(tmp_abs_ll, PREC);
    mpc_set_prec(tmp_abs_lu, PREC);
    mpc_set_prec(tmp_abs_ul, PREC);
    mpc_set_prec(tmp_abs_uu, PREC);
    
    mpc_set_r_r(tmp_abs_n, complex->re->n, complex->im->n);
    mpc_abs(tmp_abs_n_re, tmp_abs_n, MPC_RNDNN);
    
    mpc_set_r_r(tmp_abs_ll, complex->re->l, complex->im->l);
    mpc_set_r_r(tmp_abs_lu, complex->re->l, complex->im->u);
    mpc_set_r_r(tmp_abs_ul, complex->re->u, complex->im->l);
    mpc_set_r_r(tmp_abs_uu, complex->re->u, complex->im->u);
    
    int ll_mpc_ternary = mpc_abs(tmp_abs_ll_re, tmp_abs_ll, MPC_RNDDD);
    int lu_mpc_ternary = mpc_abs(tmp_abs_lu_re, tmp_abs_lu, MPC_RNDDU);
    int ul_mpc_ternary = mpc_abs(tmp_abs_ul_re, tmp_abs_ul, MPC_RNDUD);
    int uu_mpc_ternary = mpc_abs(tmp_abs_uu_re, tmp_abs_uu, MPC_RNDUU);
    
    Bounds* bounds = bounds_new();
    bounds_update_re(bounds, tmp_abs_ll_re, complex->re->L, MPC_INEX_RE(ll_mpc_ternary) == 0);
    bounds_update_re(bounds, tmp_abs_lu_re, complex->re->L, MPC_INEX_RE(lu_mpc_ternary) == 0);
    bounds_update_re(bounds, tmp_abs_ul_re, complex->re->R, MPC_INEX_RE(ul_mpc_ternary) == 0);
    bounds_update_re(bounds, tmp_abs_uu_re, complex->re->R, MPC_INEX_RE(uu_mpc_ternary) == 0);
    
    if (interval_span_zero_p(complex->im))
    {
        l_abs = real_abs(complex->re->l);
        u_abs = real_abs(complex->re->u);
        
        mpfr_set_r(tmp_abs_re_abs, real_lt_p(l_abs, u_abs) ? l_abs : u_abs);
        bounds_update_re(bounds, tmp_abs_re_abs, true, true);
        
        real_free(l_abs);
        real_free(u_abs);
    }
    
    if (interval_span_zero_p(complex->re))
    {
        l_abs = real_abs(complex->im->l);
        u_abs = real_abs(complex->im->u);
        
        mpfr_set_r(tmp_abs_im_abs, real_lt_p(l_abs, u_abs) ? l_abs : u_abs);
        bounds_update_re(bounds, tmp_abs_im_abs, true, true);
        
        real_free(l_abs);
        real_free(u_abs);
    }
    
    if (interval_span_zero_p(complex->re) && interval_span_zero_p(complex->im))
    {
        mpfr_set_zero(tmp_abs_zero, 1);
        bounds_update_re(bounds, tmp_abs_zero, true, true);
    }
    
    Interval* re = interval_new(real_from_mpfr(bounds->re_min, ROUND_DOWN), bounds->re_min_closed, real_from_mpfr(tmp_abs_n_re, ROUND_NEAREST), bounds->re_max_closed, real_from_mpfr(bounds->re_max, bounds->re_max_exact ? ROUND_UP : ROUND_UP1));
    Interval* im = interval_zero();
    
    bounds_free(bounds);
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing the area of the
 *  region in the complex plane enclosed by complex.
 */
Complex* complex_area (Complex* complex)
{
    Real* re_width = real_subtract(complex->re->u, complex->re->l);
    Real* im_width = real_subtract(complex->im->u, complex->im->l);
    Real* area = real_multiply(re_width, im_width);
    
    real_free(re_width);
    real_free(im_width);
    
    return complex_from_real(area);
}

/*
 *  Returns a pointer to a newly allocated Complex containing the length of the
 *  diagonal of the region in the complex plane enclosed by complex.
 */
Complex* complex_diagonal (Complex* complex)
{
    Real* re_width = real_subtract(complex->re->u, complex->re->l);
    Real* im_width = real_subtract(complex->im->u, complex->im->l);
    
    Real* re_width2 = real_multiply(re_width, re_width);
    Real* im_width2 = real_multiply(im_width, im_width);
    
    Real* diag2 = real_add(re_width2, im_width2);
    
    mpfr_set_r(tmp_diagonal, diag2);
    int mpfr_ternary = mpfr_sqrt(tmp_diagonal, tmp_diagonal, MPFR_RNDN);
    
    Interval* re = interval_new(real_from_mpfr(tmp_diagonal, ROUND_DOWN), true, real_from_mpfr(tmp_diagonal, ROUND_NEAREST), true, real_from_mpfr(tmp_diagonal, (mpfr_ternary == 0) ? ROUND_UP : ROUND_UP1));
    
    real_free(re_width);
    real_free(im_width);
    real_free(re_width2);
    real_free(im_width2);
    real_free(diag2);
    
    return complex_new(re, interval_zero());
}

/*
 *  Returns a pointer to a newly allocated Complex containing the conjugate of
 *  complex.
 */
Complex* complex_conjugate (Complex* complex)
{
    Interval* re = interval_dup(complex->re);
    Interval* im = interval_negate(complex->im);
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing -complex.
 */
Complex* complex_negate (Complex* complex)
{
    Interval* re = interval_negate(complex->re);
    Interval* im = interval_negate(complex->im);
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing the reflection of
 *  complex over the line re == im.
 */
Complex* complex_reflect (Complex* complex)
{
    Interval* re = interval_dup(complex->im);
    Interval* im = interval_dup(complex->re);
    
    return complex_new(re, im);
}

/*
 *  Returns a pointer to a newly allocated Complex containing √complex.
 */
Complex* complex_sqrt (Complex* complex)
{
    mpfr_set_prec(tmp_sqrt_n_re, PREC);
    mpfr_set_prec(tmp_sqrt_n_im, PREC);
    mpc_set_prec(tmp_sqrt_n, PREC);
    mpc_set_prec(tmp_sqrt_ll, PREC);
    mpc_set_prec(tmp_sqrt_lu, PREC);
    mpc_set_prec(tmp_sqrt_ul, PREC);
    mpc_set_prec(tmp_sqrt_uu, PREC);
    mpc_set_prec(tmp_sqrt_re_sqrt, PREC);
    
    mpc_set_r_r(tmp_sqrt_n, complex->re->n, complex->im->n);
    mpc_sqrt(tmp_sqrt_n, tmp_sqrt_n, MPC_RNDNN);
    mpc_real(tmp_sqrt_n_re, tmp_sqrt_n, MPFR_RNDN);
    mpc_imag(tmp_sqrt_n_im, tmp_sqrt_n, MPFR_RNDN);
    
    mpc_set_r_r(tmp_sqrt_ll, complex->re->l, complex->im->l);
    mpc_set_r_r(tmp_sqrt_lu, complex->re->l, complex->im->u);
    mpc_set_r_r(tmp_sqrt_ul, complex->re->u, complex->im->l);
    mpc_set_r_r(tmp_sqrt_uu, complex->re->u, complex->im->u);
    
    int ll_mpc_ternary = mpc_sqrt(tmp_sqrt_ll, tmp_sqrt_ll, MPC_RNDDD);
    int lu_mpc_ternary = mpc_sqrt(tmp_sqrt_lu, tmp_sqrt_lu, MPC_RNDDU);
    int ul_mpc_ternary = mpc_sqrt(tmp_sqrt_ul, tmp_sqrt_ul, MPC_RNDUD);
    int uu_mpc_ternary = mpc_sqrt(tmp_sqrt_uu, tmp_sqrt_uu, MPC_RNDUU);
    
    Bounds* bounds = bounds_new();
    bounds_update(bounds, tmp_sqrt_ll, complex->re->L, complex->im->L, ll_mpc_ternary);
    bounds_update(bounds, tmp_sqrt_lu, complex->re->L, complex->im->R, lu_mpc_ternary);
    bounds_update(bounds, tmp_sqrt_ul, complex->re->R, complex->im->L, ul_mpc_ternary);
    bounds_update(bounds, tmp_sqrt_uu, complex->re->R, complex->im->R, uu_mpc_ternary);
    
    if (interval_span_zero_p(complex->im))
    {
        mpc_set_r_r(tmp_sqrt_re_sqrt, complex->re->l, 0);
        int re_sqrt_mpc_ternary = mpc_sqrt(tmp_sqrt_re_sqrt, tmp_sqrt_re_sqrt, MPC_RNDDD);
        
        bounds_update(bounds, tmp_sqrt_re_sqrt, complex->re->L, true, re_sqrt_mpc_ternary);
    }
    
    Interval* re = interval_new(real_from_mpfr(bounds->re_min, ROUND_DOWN), bounds->re_min_closed, real_from_mpfr(tmp_sqrt_n_re, ROUND_NEAREST), bounds->re_max_closed, real_from_mpfr(bounds->re_max, bounds->re_max_exact ? ROUND_UP : ROUND_UP1));
    Interval* im = interval_new(real_from_mpfr(bounds->im_min, ROUND_DOWN), bounds->im_min_closed, real_from_mpfr(tmp_sqrt_n_im, ROUND_NEAREST), bounds->im_max_closed, real_from_mpfr(bounds->im_max, bounds->im_max_exact ? ROUND_UP : ROUND_UP1));
    
    bounds_free(bounds);
    
    return complex_new(re, im);
}
